%% Risk Neutral Marke Maker
clc;
clear;

N_M=2;
N_I=10;


lambda_I =1;


T=20;
tau_V=1;

tau_eps=1*ones(T,1);


tau_eta=10*ones(T,1);




theta_bar=1;
V_bar=1;   

X_0=0;


theta_I_bar=theta_bar/N_I;




o_I_V=zeros(T,1);
for j=1:T
    o_I_V(j)=1/(tau_V+sum(tau_eps(1:j)));
end
K_I_V=o_I_V.*tau_eps;

sigma_I=zeros(T,1);     
sigma_I(1)=1/tau_V+1/tau_eps(1); 
for j=2:T
    sigma_I(j)=o_I_V(j-1)+1/tau_eps(j);  
end

Sigma_I=zeros(2,2,T);
for t=1:T
    Sigma_I(1,1,t)=sigma_I(t);
    Sigma_I(2,2,t)=1/tau_eta(t);
end



% mu=zeros(T,1);
% mu(T)=lambda_I*o_I_V(T);
% for t=T:-1:2
%     mu(t-1)=mu(t)+lambda_I*sigma_I(t)*(K_I_V(t))^2;
% end

mu= lambda_I*o_I_V;


h=lambda_I./tau_eps;

Sigma_U=zeros(T,1);
Sigma_U(1)=1/tau_V+1/tau_eps(1)+h(1)^2/tau_eta(1);

o_U_V=zeros(T,1);
o_U_X=zeros(T,1);
o_U_VX=zeros(T,1);

o_U_V(1)=(tau_eps(1)^-1+h(1)^2*tau_eta(1)^-1)/(tau_V*Sigma_U(1));
o_U_X(1)=1/tau_eta(1)-h(1)^2/(tau_eta(1)^2*Sigma_U(1));
o_U_VX(1)=(1/tau_V)*(h(1)/tau_eta(1))/Sigma_U(1);

K_U_V=zeros(T,1);
K_U_X=zeros(T,1);

K_U_V(1)=tau_V^-1/Sigma_U(1);
K_U_X(1)=-h(1)*tau_eta(1)^-1/Sigma_U(1);

for j=2:T
    Sigma_U(j)=o_U_V(j-1)+1/tau_eps(j)+h(j)^2/tau_eta(j);
    
    o_U_V(j)=(1/tau_eps(j)+h(j)^2/tau_eta(j))*o_U_V(j-1)/Sigma_U(j);
    o_U_X(j)=o_U_X(j-1)+1/tau_eta(j)-(o_U_VX(j-1)-h(j)/tau_eta(j))^2/Sigma_U(j);
    o_U_VX(j)=o_U_VX(j-1)-(o_U_VX(j-1)-h(j)/tau_eta(j))*o_U_V(j-1)/Sigma_U(j);
    
    K_U_V(j)=o_U_V(j-1)/Sigma_U(j);
    K_U_X(j)=(o_U_VX(j-1)-h(j)/tau_eta(j))/Sigma_U(j);
end

H_I=zeros(4,4,T);
for t=2:T
    H_I(1,1,t)=1;
    H_I(2,2,t)=1;
    H_I(3,3,t)=1-K_U_X(t)*mu(t-1);
    H_I(3,2,t)=K_U_X(t)*mu(t-1);
    H_I(4,4,t)=1;
end

F_I=zeros(4,2,T);
for t=1:T
    F_I(1,1,t)=K_I_V(t);
    F_I(2,2,t)=1;
    F_I(3,1,t)=K_U_X(t);
    F_I(3,2,t)=-K_U_X(t)*h(t);
end

F_M=zeros(3,1,T);
for t=1:T
    F_M(1,1,t)=K_U_V(t);
    F_M(2,1,t)=K_U_X(t);
end








% the values of the parameters for the last period

gamma_I=zeros(T,1);
gamma_I(T)=lambda_I*o_I_V(T);


delta_I=zeros(T,1);
delta_I(T)=1;

delta_M=zeros(T,1);
delta_M(T)=0;

g_I=zeros(T,1);
g_I(T)=N_M*lambda_I*o_I_V(T)/(N_M+1);

g_M=zeros(T,1);
g_M(T)=-N_I/(N_M+1);

f_I=zeros(T,1);
f_I(T)=lambda_I*o_I_V(T)/(N_M+1);

f_M=zeros(T,1);
f_M(T)=N_I/(N_M+1);

rho_I=zeros(T,1);
rho_I(T)=1;

rho_M=zeros(T,1);
rho_M(T)=0;

m_I=zeros(T,1);
m_I(T)=lambda_I*o_I_V(T)/(N_M+1)^2;

m_M=zeros(T,1);
m_M(T)=2*N_I*lambda_I*o_I_V(T)/(N_M+1)^2;




C_I=zeros(4,1,T);
C_I(3,1,T)=-lambda_I*o_I_V(T)*N_M/(N_M+1)^2;


C_M=zeros(3,1,T);
C_M(2,1,T)=2*N_I*lambda_I*o_I_V(T)/(N_M+1)^2;

M_I=zeros(4,4,T);
M_I(1,2,T)=1;
M_I(2,1,T)=1;
M_I(2,2,T)=-lambda_I*o_I_V(T);
M_I(2,4,T)=-1;
M_I(3,3,T)=lambda_I*o_I_V(T)*(N_M/(N_M+1))^2;
M_I(4,2,T)=-1;

M_M=zeros(3,3,T);
M_M(2,2,T)=2*N_I*lambda_I*o_I_V(T)/(N_M+1)^2;



Xi_I=zeros(2,2,T);
Xi_I(:,:,T)=eye(2)/( eye(2)/Sigma_I(:,:,T)+ lambda_I*transpose(F_I(:,:,T))*M_I(:,:,T)*F_I(:,:,T));
% Xi_I(:,:,T)=inv(inv(Sigma_I(:,:,T))+ lambda_I*transpose(F_I(:,:,T))*M_I(:,:,T)*F_I(:,:,T));



H_I_P=zeros(4,1,T);
H_I_P(1,1,T)=delta_I(T);
H_I_P(2,1,T)=-mu(T)*delta_I(T);
H_I_P(3,1,T)=g_I(T);
H_I_P(4,1,T)=1-delta_I(T);

H_M_alpha=zeros(3,1,T);
H_M_alpha(1,1,T)= -delta_M(T);
H_M_alpha(2,1,T)= mu(T)*delta_M(T)+g_M(T);
H_M_alpha(3,1,T)= delta_M(T);

H_I_alpha=zeros(4,1,T);
H_I_alpha(1,1,T)= -delta_M(T);
H_I_alpha(2,1,T)= delta_M(T)*mu(T);
H_I_alpha(3,1,T)= g_M(T);
H_I_alpha(4,1,T)= delta_M(T);

H_I_theta=zeros(4,1,T);
H_I_theta(:,:,T)=H_I_alpha(:,:,T)*N_M/N_I;

f_I_theta=zeros(T,1);
f_I_theta(T)=1-f_M(T)*N_M/N_I;




H_I_D=zeros(4,1,T);







test1=zeros(T,1);


test2=zeros(T,1);
test2(T)=1;

test3=zeros(T,1);


SOD_I=zeros(T-1,1);
SOD_M=zeros(T-1,1);

q_M=zeros(T-1,1);

omega=zeros(T,1);
omega(T)=N_M/(N_M+1);

% coefficients in P_MR
gamma_MR=zeros(T,1);  
delta_MR=zeros(T,1);
delta_MR(T)=1;
g_MR=zeros(T,1);

% coefficients in P_IR
delta_IR=zeros(T,1);
delta_IR(T)=1;
g_IR=zeros(T,1);


% compute the values of the parameters recursively

for t=T:-1:2
    
    
    H_I_D(:,:,t)=transpose( transpose(H_I_P(:,:,t))*H_I(:,:,t)-...
        lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)...
        *transpose(F_I(:,:,t))*M_I(:,:,t)*H_I(:,:,t) );

    delta_IR(t-1)=H_I_D(1,1,t);
    g_IR(t-1)=H_I_D(3,1,t);
    
    
    test1(t-1)=-H_I_D(2,1,t)/H_I_D(1,1,t)-mu(t-1);   % =mu(t-1)
    test2(t-1)=H_I_D(1,1,t)+H_I_D(4,1,t);  % =1
    
       
    gamma_I(t-1)=f_I(t)+lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)...
        *transpose(F_I(:,:,t))*(H_I_P(:,:,t)+C_I(:,:,t));
    
    SOD_I(t-1)=lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*H_I_P(:,:,t);
    SOD_M(t-1)=gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2;

    q_M(t-1)=SOD_M(t-1);
    omega(t-1)=N_M*gamma_I(t-1)/(N_I*q_M(t-1));

    gamma_MR(t-1)=m_M(t)/N_I;
    delta_MR(t-1)=1-C_M(1,1,t)/N_I;
    g_MR(t-1)=C_M(2,1,t)/N_I;
    
    
    delta_M(t-1)=( 1-H_I_D(1,1,t)-C_M(1,1,t)/N_I )/ (gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2);
    
    g_M(t-1)=( H_I_D(3,1,t)+C_M(2,1,t)/N_I+ mu(t-1)*C_M(1,1,t)/N_I-mu(t-1) )...
        /(gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2);
    
    f_M(t-1)=( gamma_I(t-1)-m_M(t)/N_I )/(gamma_I(t-1)*(N_M+1)/N_I-m_M(t)*N_M/(N_I)^2);
   
    f_I_theta(t-1)=1-f_M(t-1)*N_M/N_I;
    
    H_I_alpha(1,1,t-1)= -delta_M(t-1);
    H_I_alpha(2,1,t-1)= delta_M(t-1)*mu(t-1);
    H_I_alpha(3,1,t-1)= g_M(t-1);
    H_I_alpha(4,1,t-1)= delta_M(t-1);
    
    H_M_alpha(1,1,t-1)=-delta_M(t-1);
    H_M_alpha(2,1,t-1)=delta_M(t-1)*mu(t-1)+g_M(t-1);
    H_M_alpha(3,1,t-1)=delta_M(t-1);

    H_I_theta(:,:,t-1)=H_I_alpha(:,:,t-1)*N_M/N_I;
    
    delta_I(t-1)=H_I_D(1,1,t)+gamma_I(t-1)*delta_M(t-1)*N_M/N_I;
    
    g_I(t-1)=H_I_D(3,1,t)-gamma_I(t-1)*g_M(t-1)*N_M/N_I;
    
    f_I(t-1)=gamma_I(t-1)-gamma_I(t-1)*f_M(t-1)*N_M/N_I;
    
    H_I_P(1,1,t-1)=delta_I(t-1);
    H_I_P(2,1,t-1)=-mu(t-1)*delta_I(t-1);
    H_I_P(3,1,t-1)=g_I(t-1);
    H_I_P(4,1,t-1)=1-delta_I(t-1);
    
    
    
    
    
    m_I(t-1)= (f_I_theta(t-1))^2* (m_I(t)+lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t)- lambda_I*transpose(C_I(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*C_I(:,:,t) );
      
    
    m_M(t-1)=m_M(t)*(1-2*f_M(t-1)*(N_M-1)/N_I) +...
        f_M(t-1)^2* (2*gamma_I(t-1)/N_I+ m_M(t)* (N_M-2)*N_M/(N_I)^2);

  
    
   
    
    C_I(:,:,t-1)=f_I_theta(t-1)*( transpose(H_I(:,:,t))*C_I(:,:,t)+ m_I(t)*H_I_theta(:,:,t-1)+...
        lambda_I* H_I_theta(:,:,t-1)*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*H_I_P(:,:,t)- lambda_I*f_I_theta(t-1)*transpose( C_I(:,:,t)*...
        transpose(H_I_theta(:,:,t-1))+ transpose(M_I(:,:,t))*H_I(:,:,t) )*F_I(:,:,t)*Xi_I(:,:,t)*...
        transpose(F_I(:,:,t))*C_I(:,:,t) );
    
    C_M(:,:,t-1)= (1-f_M(t-1)*(N_M-1)/N_I)*C_M(:,:,t)+...
        H_M_alpha(:,:,t-1)*...
        (m_M(t)*(N_M-1)/N_I- m_M(t)*f_M(t-1)*(N_M-2)*N_M/(N_I)^2-2*f_M(t-1)*gamma_I(t-1)/N_I );
    
    test3(t-1)=C_M(1,1,t-1)+C_M(3,1,t-1); %=0
    
    
    
    M_I(:,:,t-1)=transpose(H_I(:,:,t))*M_I(:,:,t)*H_I(:,:,t)+m_I(t)*H_I_theta(:,:,t-1)*...
        transpose(H_I_theta(:,:,t-1))+ H_I_theta(:,:,t-1)*transpose(C_I(:,:,t))*H_I(:,:,t)+...
        transpose(H_I_theta(:,:,t-1)*transpose(C_I(:,:,t))*H_I(:,:,t))+...
        lambda_I*transpose(H_I_P(:,:,t))*F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*...
        H_I_P(:,:,t)*H_I_theta(:,:,t-1)*transpose(H_I_theta(:,:,t-1))-...
        lambda_I*transpose( C_I(:,:,t)*transpose(H_I_theta(:,:,t-1))+ transpose(M_I(:,:,t))*H_I(:,:,t) )*...
        F_I(:,:,t)*Xi_I(:,:,t)*transpose(F_I(:,:,t))*( C_I(:,:,t)*...
        transpose(H_I_theta(:,:,t-1))+ transpose(M_I(:,:,t))*H_I(:,:,t) );
    
    M_M(:,:,t-1)=M_M(:,:,t)+ (2*gamma_I(t-1)/N_I+ m_M(t)*(N_M-2)*N_M/(N_I)^2)*H_M_alpha(:,:,t-1)*...
        transpose(H_M_alpha(:,:,t-1)) + H_M_alpha(:,:,t-1)*transpose(C_M(:,:,t))*(N_M-1)/N_I+...
        C_M(:,:,t)*transpose(H_M_alpha(:,:,t-1))*(N_M-1)/N_I;
    
    
    
    Xi_I(:,:,t-1)= eye(2)/(eye(2)/Sigma_I(:,:,t-1)+...
        lambda_I*transpose(F_I(:,:,t-1))*M_I(:,:,t-1)*F_I(:,:,t-1));

    
    rho_I(t-1)=rho_I(t)*sqrt(det(Xi_I(:,:,t))/det(Sigma_I(:,:,t)) );
    
    rho_M(t-1)=rho_M(t)+ transpose(F_M(:,:,t))*M_M(:,:,t)*F_M(:,:,t)*Sigma_U(t)/2;   

end

rho_M_0=rho_M(1)+ transpose(F_M(:,:,1))*M_M(:,:,1)*F_M(:,:,1)*Sigma_U(t)/2; 

SOC_I=find(SOD_I<0)
SOC_M=find(SOD_M<0)








figure(1)
plot(0:1/(T-1):1, delta_IR-(g_IR+g_MR)./mu, 0:1/(T-1):1, delta_MR-(g_IR+g_MR)./mu, '-.')
xlabel('$t/T$','interpreter','latex')
legend('$\xi^I_t-\xi^U_t$', '$\psi^U_t-\psi^I_t$', 'interpreter','latex')



